Noname manuscript No. 

(will be inserted by the editor) 



DEBRIS DISK SIZE DISTRIBUTIONS: STEADY STATE 
COLLISIONAL EVOLUTION WITH P R DRAG AND 
OTHER LOSS PROCESSES 

M. C. WYATT 1 , C. J. CLARKE 1 and M. BOOTH 1 2 3 

1 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 OHA, 

UK 

E-mail: wyatt@ast.cam.ac.uk 
2 Herzberg Institute of Astrophysics, 5071 West Saanich Road, Victoria, B. C, V9E 2E7, 

Canada 

■ 3 University of Victoria, Finnerty Road, Victoria, B. C, V8W 3P6, Canada 



oo 
(N 



O 



Submitted: 26 January 2011, Revised: 1 March 2011 



Abstract We present a new scheme for determining the shape of the size distribution, and 
its evolution, for collisional cascades of planetesimals undergoing destructive collisions and 
loss processes like Poynting-Robertson drag. The scheme treats the steady state portion of 
the cascade by equating mass loss and gain in each size bin; the smallest particles are ex- 
Oh ■ pected to reach steady state on their collision timescale, while larger particles retain their 

primordial distribution. For collision-dominated disks, steady state means that mass loss 
rates in logarithmic size bins are independent of size. This prescription reproduces the ex- 
pected two phase size distribution, with ripples above the blow-out size, and above the tran- 
sition to gravity-dominated planetesimal strength. The scheme also reproduces the expected 
evolution of disk mass, and of dust mass, but is computationally much faster than evolving 
distributions forward in time. For low-mass disks, P-R drag causes a turnover at small sizes 
to a size distribution that is set by the redistribution function (the mass distribution of frag- 
ments produced in collisions). Thus information about the redistribution function may be 
recovered by measuring the size distribution of particles undergoing loss by P-R drag, such 
as that traced by particles accreted onto Earth. Although cross-sectional area drops with age 
Q\ . °c t~ 2 in the PR-dominated regime, dust mass falls <x f~ 2 - 8 , underlining the importance of 

' understanding which particle sizes contribute to an observation when considering how disk 

detectability evolves. Other loss processes are readily incorporated; we also discuss gener- 
alised power law loss rates, dynamical depletion, realistic radiation forces and stellar wind 
CT) | drag. 



Keywords circumstellar matter • planetary systems • debris disks • collisions 



1 Introduction 

X 

Many nearby main sequence stars are seen to be surrounded by /xm-sized dust known as 
debris disks (see review in Wyatt 2008). The dust is short-lived indicating that larger plan- 
etesimals are present that feed the observed dust (Backman & Paresce 1993). It is thought 
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that the large planetesimals are destroyed in mutual collisions, creating fragments that also 
collide to produce yet smaller fragments, and so on until the dust is small enough to be re- 
moved from the system by radiation pressure; i.e., that there is a continuous size distribution 
extending from small to large objects known as a collisional cascade. 

In this respect extrasolar debris disks are similar to the asteroid and Kuiper belts in the 
Solar System. The size distribution of large objects in the Solar System's belts is well char- 
acterised observationally (Gladman et al. 2009; Fuentes et al. 2010). There remains debate 
as to the origin of this distribution at the largest sizes (Durda et al. 1998; Bottke et al. 2005; 
Morbidelli et al. 2009), but it is recognised that these populations have undergone, and con- 
tinue to undergo, collisional evolution which means that their size distributions should also 
extend down to dust sizes. Exactly how the distributions extrapolate down is less well con- 
strained (e.g., Schlichting et al. 2009). However, dust is seen in the inner Solar System that 
migrates inwards past the Earth due to Poynting-Robertson drag (Leinert & Griin 1990; Der- 
mott et al. 2001), though it is debated as to whether this dust has an asteroidal or cometary 
origin (Durda & Dermott 1997; Nesvorny et al. 2010). 

One big difference with extrasolar debris disks is that, to be detected above the light 
from the star, their dust has to be orders of magnitude more abundant than dust in the Solar 
System. This means that it has been possible to ignore the effects of P-R drag when studing 
their evolution, since for such disks collisional timescales (that scale inversely with total 
dust mass) are much shorter than P-R drag timescales (that are independent of dust mass) 
(Wyatt 2005). The regime of such collision-dominated disks is now (relatively) well under- 
stood. Early analytical models treated systems in which the size distribution is defined by a 
single power law, and showed that disk luminosity is expected to stay constant, but to fall off 
inversely with age at late times (Dominik & Decin 2003; Wyatt et al. 2007a), in a manner 
consistent with the observed debris disk population (Wyatt et al. 2007b). More recent analyt- 
ical models allow for three phases in the size distribution to allow for planetesimals having 
a size-dependent dispersal threshold (Lohne et al. 2008; Heng & Tremaine 2009), and pro- 
vide a better fit to the size distribution expected from more detailed numerical models that 
solve for the evolution of an input size distribution given assumptions about the outcome of 
collisions between different sized particles (Krivov et al. 2006; Thebault & Augereau 2007). 
However, the analytical models cannot fit details in the expected size distribution, which is 
not just defined by power laws, and neither do they consider effects other than collisions. 

There is now a need to understand how loss processes affect the size distribution, and 
so the detectability, of a debris disk. Current instrumentation is probing dust levels at which 
drag may start to become important, and as searches for Earth-analogs intensify, it is in- 
evitable that lower levels of dust will be probed and the dynamics in this regime needs to 
be properly understood (Beichman et al. 2006). It is also becoming clear that stellar wind 
drag, which acts in a similar manner to P-R drag, may be important for young stars with 
high mass loss rates (Plavchan et al. 2005; Reidemeister et al. 2010). 

Although the competition between collisions and P-R drag has been considered in the 
Solar System (Leinert et al. 1983; Griin et al. 1985; Ishimoto 2000), such studies focus on 
explaining the currently observed dust distribution, and are not readily applicable to earlier 
epochs or to the debris disks of other stars. Collisional models that do study the long-term 
evolution of the asteroid belt typically exclude small particles and focus on fitting the large 
body size distribution (Bottke et al. 2005), whereas models for structures in the zodiacal 
cloud often only treat the dynamics of individual particles essentially ignoring collisions 
(Dermott et al. 1994; Grogan et al. 2001). However, more recently the balance of collisions 
and P-R drag has been incorporated into codes that evolve size distributions with time and 
so far results have been presented for two systems (Vitense et al. 2010; Reidemeister et al. 
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2010). Other models are also being developed that incorporate a (less detailed) prescription 
for collisions into N-body simulations (Stark & Kuchner 2009; Kuchner & Stark 2010). 

Here we present a new scheme that can be used to determine the size distribution of a 
planetesimal belt undergoing destructive collisions as well as loss processes. This scheme is 
not expected to provide a more faithful measure of the size distribution than would be found 
by evolving a distribution forward in time. Rather its value is in the speed with which the 
distribution can be determined (in seconds), with its ability to reproduce details in the steady 
state size distribution, including when loss processes are acting, and with the simplicity of 
the analytical presentation which permits an analytical understanding of how various pro- 
cesses affect the size distribution. Although the scheme presented here is in its most basic 
form, considering the 1-dimensional size distribution (rather than, e.g., including spatial di- 
mensions), and treating collisional destruction with a 1-dimensional redistribution function 
(essentially averaging the outcome over all possible collisions), extension to extra dimen- 
sions is possible in principle. |J2] outlines the development of the scheme by application to 
the collision-dominated regime, then jj3] continues by detailing the effect of P-R drag; |j4] 
considers other loss processes, and conclusions are given in Jj5] 



2 Steady state in the collisional regime 

Consider a belt of planetesimals in which is the total mass in the A:-th bin, where is the 
typical size of planetesimals in the bin, and bins are logarithmically spaced in size, with k = 1 
being the largest bin and working to smaller sizes with increasing k so that D^i/D^ = 1 — 8 
where 8 is approximately the logarithmic bin width, since it is also assumed that 8 <C 1. 
Collisions between planetesimals cause the mass to be redistributed amongst the bins, and 
furthermore loss processes may cause mass to be removed from the belt, and there may be 
sources that add mass to the belt. 

From mass conservation for the k-th bin 

m k = m k -m k -m k + m k s , (1) 

where the superscripts +c, — c, —I, +g refer to the mass gained from collisions in other bins, 
the mass lost from collisions in this bin, the mass removed by loss processes, and the mass 
gained from other sources, respectively. Note that our 1-dimensional prescription of the belt 
does not have a radial dimension and so considers loss processes such as P-R drag as a 
sink, rather than as a diffusive process (e.g. Gorkavyi et al. 1997). For now we consider the 
situation where rht 8 = m k l = 0, but return to loss processes in [J3] 



2. 1 Mass loss is independent of size 

We define the redistribution function f(i,k) to be the fraction of mass leaving the i-th bin 
from collisions that goes into the k-th bin. Since mass must be conserved JjjLj f(i, k) = 1. In 
this paper the model is being applied to situations in which collisions are destructive, so that 
f(i, k < i) = and f(hfy = 1» though it may be possible to extend this to situations 

where accretion can also take place (see, e.g., Birnstiel et al. 2011). 

Let us further assume that the redistribution function is scale independent, so that /(«' + 
n,k + n) = f(i,k), where n is any integer, which we can implement by defining a new redis- 
tribution function needing just one parameter F(k—i)=f(i,k). Putting this into equation (TJ 
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and giving collisional gains and losses in terms of the redistribution function we find 

k-\ 

m k =£m; c F(k-i)-m^. (2) 

If we assume that the size distribution is in quasi-steady state, so that rh^ w 0, and append 
s to any subscripts of quantities that are calculated in steady state, we find that 

*b=E^(*-0- (3) 

Since YT=i ^(0 = ^> men as l° n § as TT=k^0) ^> one solution to equation <[3j> must be 
that 

'V = '«,r = C, (4) 

where C is a constant; i.e., the mass loss rate in logarithmic bins is independent of size. 

Note that the main assumption to arrive at this result was that the redistribution function 
is scale independent. However, there is a further assumption regarding whether the size 
distribution is infinite in extent (e.g., Dohnanyi 1969). Truncated distributions are considered 
in more detail in following sections, but here we note that for distributions truncated at 
large sizes, equation (|4j is expected to apply at sizes for which the mass that would have 
been replenished from sizes above the truncation is negligible (i.e., it applies for k 2> 1 and 
redistribution functions that are weighted toward larger sizes). 



2.2 Mass loss rate 

Equations lO and © can be used to determine the steady state size distribution, m^. To do 
so, we consider that the mass loss rate from the k-th bin due to collisions is 

m- c = mk R c k) (5) 

where R c k is the rate at which individual objects of size undergo catastrophic collisions. 

A catastrophic collision is defined as one in which the mass of the largest remnant after 
the collision has less than half of the mass of the original object. This definition is used 
for our typical collision, not because we believe that catastrophic collisions are the only 
important collisions — indeed it has been shown that cratering collisions with much lower 
energy may also result in significant mass loss (e.g., Kobayashi & Tanaka 2010) — but 
because we expect the overall mass loss rate (i.e., including both catastrophic and cratering 
collisions) for all bins to scale with the rate of catastrophic collisions within them, at least 
for distributions that are infinite in extent. It may be that mass loss in fact scales with the 
rate of collisions in which the largest remnant has slightly more or less than 50% the mass 
of the original object. 

The catastrophic collision rate for a planetesimal of size can be determined using the 
particle-in-a-box approach and integrating over the size distribution of impactors that have 
sufficient specific incident energy (Q) to cause a catastrophic collision, which is defined as 
the dispersal threshold Q* D . For collisions that occur at relative velocity v re / the smallest 
impactors that cause catastrophic destruction are those of size X c Dk, where 



(6) 
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which we denote in the following equation as those with an index i ck . Thus the collision rate, 
both in its discrete and continuous forms, is 



3wj 



^E^^+AO^, (7) 



f Di n{D i )Q.25{D k +D i fPi k dD i , 

JX c D k 



(8) 



where p is particle density, P% is the intrinsic collision probability used by other authors 
(e.g., Wetherill 1967) between particles /' and k that is equal to Kv re i/V for all particle sizes 
for applications in this paper, and V is the volume through which the planetesimals are 
moving, n(D) is the number of planetesimals in the size distribution per unit diameter, and 
the gravitational focussing factor has been neglected for this calculation. 



2.3 Analytical steady state method 

The continuous form of the collision rate (equation [8} can be used to determine the steady 
state size distribution analytically, for certain assumptions. By integrating the mass distribu- 
tion (TipD 3 /6)n(D) over the bin we find that for small 8 

m k « (K P /6)n(D k )D 4 k 8, (9) 

and so 

m^ocn(D k )D 4 k [ ' n{D t )(D k + D^dD,. (10) 

Jx c D k 

Further assuming that 



n(D)=KD- a , (11) 
and that a > 3 and X c <C 1 so that collisions close to the catastrophic collision threshold 



dominate the rate (i.e., so we can take only the term with D 2 k in the integral and ignore the 



upper limit), gives 
HQ* D <*D b then we find that 



x l-a D l-2a (12) 



m k <*D k , (13) 

and for this to be independent of size as required for steady state (equation [4} 

a = (7 + 6/3)/(2 + fe/3). (14) 

This reduces to the canonical a = 3.5 distribution for the situation where strength is inde- 
pendent of size (Dohnanyi 1969; Tanaka et al. 1996), and replicates the distribution found 
analytically by other authors when strength depends on size (O'Brien & Greenberg 2003; 
Lohne 2008; Lohne et al. 2008; Heng & Tremaine 2009; Birnstiel et al. 2011; Belyaev & 
Rafikov2011). 



2.4 Numerical steady state method 



The discrete form of the collision rate (equation |7} can also be used to determine the steady 
state size distribution numerically. 
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2.4.1 Method 1 

If we impose a constant mass loss rate per logarithmic bin, as argued to be the case in £12.11 
then we can rearrange equation (|5} and combine with equation I0 to get 

m ks = C/R c k . (15) 

This can be solved iteratively, and converges quickly to find the steady state distribution. 
The shape of the resulting distribution is independent of the total mass in the distribution, 
because multiplying the whole distribution by TV would result in the mass loss in the 
k-th bin m^ c and the mass gain rh^ c both increasing by TV 2 . Application of this method is 
discussed in £12.6.21 

2.4.2 Method 2 

Alternatively we can relax the assumption of equation © and apply a similar method more 
directly to equation ifJJ, and so equate mass loss with mass gain in each bin by iterating 

m ks = C k /R c k , (16) 

k-1 

Ck=Y.' i hs C F(k-i). (17) 

i=i 

In this case, by the arguments of £12.11 we would expect to end up with a solution for which 
Ck is approximately constant (at least for most particle sizes). 

To implement this method requires the redistribution function F{1) to be defined (see 
£]2.5t . Without further constraints this method would converge to a distribution with zero 
mass. This is because the lack of larger particles means that the mass input rate to the top 
bin is zero, which can only be balanced by having zero mass there; if that is the case then 
the mass input to the bin below this is zero, and so on. This issue can be remedied by con- 
straining the largest bin(s) to have either a fixed mass distribution, or a fixed mass input rate 
(e.g., by adding rii^ 8 to C k ). Further numerical tricks are also required to achieve conver- 
gence. If numeric subscripts indicate the iteration number, or more specifically a variable 
calculated using the mass distribution from that iteration, then we found convergence with 
the following two-step scheme 

mfai.5 = [(C k i/R c kl )+m ksl }/2, (18) 
m ks 2 = [(Cki/R c kl . 5 )+mk sl }/2. (19) 

This method is employed in £12.91 and further extended to incorporate loss due to P-R drag 
in El 

2.5 Redistribution function 

Where it is necessary in this paper, we assume that the fragments from the destruction of 
objects of size Dk follow the power law of equation dllt . with parameters identified by 
subscript r, from sizes of rj rmax Dk down to infinitesimally small particles. Unless otherwise 
stated we consider r\, max = 2~ 1//3 (i.e., the largest objects have half the mass of the original), 
which corresponds to catastrophic collisions dominating the redistribution of mass. In this 
paper we also consider a function with r\ rmax = 1, which corresponds to cratering collisions 
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dominating mass redistribution. For r\ rmax < 1 we also include a small but finite fraction of 
mass in particles in the r\ rma xDk to Dk size range using a steep power law (with a = 
The slope is always assumed to have a r < 4 so that the mass in fragments is weighted 
toward large objects. These assumptions mean that the fractional mass distribution of these 
fragments is m(D) w (4 — a, )Z) 3_a ' {T\ rmax Di l ) a '^ A across the appropriate range. Integrating 
over the size range of bin D^+i, and assuming 5< 1, gives 

F(l) « <„- 4 (4- a r )S(l - 5f A - a '\ (20) 

Although it is the redistribution function F(k — i) that is employed throughout this pa- 
per, it is worth considering a more general form of the redistribution function, since the 
techniques presented herein may also applicable in this instance. Most numerical models 
of collisional evolution define the redistribution function as fz(i,j,k), the fraction of mass 
leaving the i'-th bin from collisions with the /'-th bin that goes into the k-th bin. In general 
such redistribution functions also have some form of scale independence, and can, e.g., be 
written as Fi(j — j c i,k — i), where /„■ is the index of the size that causes catastrophic de- 
struction of objects of size £),-. This is because a collision with the /' c ,-th bin ends up, by 
definition, with a largest remnant that has half the mass of the original, and collisions with 
other bins end up with largest remnants that scale as a function of the ratio Q/Q* D , and so 
have redistribution functions that only depend on j — /'„-. 

If we assume that F2 (j — j c i, k — i) is the true redistribution function, then we have im- 
plicitly averaged over j, weighted by the frequency of such collisions, to get F(k — i) as the 
redistribution of mass in a typical collision. Thus we can write 

vll L7=i>»j(Di + Dj) 2 D?F2U-j ci ,k-i) 

F(k — i) = — : . (21) 

^UmjiDt+DjfDf 

This illustrates that the F(k — i) redistribution function may not be scale independent, even 
if F2U — jci,k — i) is, because it includes information about the size distribution (see also 
Belyaev & Rafikov 201 1). In particular, it may be expected to be different at locations where 
there are sharp transitions in the size distribution. These issues may also contribute to the 
timescale of mass loss in a typical collision being different to the timescale for catastrophic 
collisions, as discussed in £12.21 Nevertheless, we retain the F(k — i) redistribution function 
for this paper, since it is expected to be valid over most of the size distribution, and is shown 
in later sections to be in agreement with more detailed models. 



2.6 Application to realistic dispersal laws 

As reported by several authors studying collision outcomes (e.g., Durda et al. 1998; Benz & 
Asphaug 1999), the size dependence of the dispersal threshold can be expressed as 

Q*D = Q a D- a + Q b D h , (22) 

where the subscripts a and b refer to contributions from the planetesimal's material strength 
and from its self-gravity, respectively, and a and b are positive. The weakest planetesimals 
then have size 

„ aN .«„+» 



°- '.*£.' ; <23) 



This is necessary to allow the size distribution to be calculated at late times in H2.9l when D, > T] rmax D\ 
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Table 1 Summary of (some of) the parameters used in the paper: power law indices that refer to equation Jilt , 
and diameters. Diameters (e.g., D x ) also have corresponding indices (e.g., i x ) to denote bin number. 



Power law index 


Regime 


Equat 


on 


a p 


Primordial 


given 


a,, 


Collisional (gravity) 


M 


a u 


Collisional (strength) 




a pr 


PR-dominated 




a, 


General loss-dominated 




53 




a r 


Redistribution function 


§ 


iven 


Diameter 


Meaning (or transition) 


Equation 


Di 


Largest planetesimals 


given 


D, 


Primordial-collisional 




29 




D w 


Strength-gravity scaling 




21 




D pr 


Collisional-PR-dominated 




37 




Di 


Collisional-loss-dominated 




5T 




D b i 


Radiation pressure-dominated 




26 





e.g., for planetesimals made of basalt impacting at 3 km s (Benz & Asphaug 1999), rea- 
sonable values are Q a = 790 J kg -1 , a = 0.38, Q b = 0.017 J kg" 1 , b = 1.36, for which 
D w = 0.23 km. Here we adopt slightly different values of Q a = 620 J kg -1 , a = 0.3, 
Qb = 5.6 x 10~ 3 J kg -1 , b = 1.5, to allow comparison with the results of Lohne et al. (2008) 
in later sections. 

2.6.1 Analytical 2-phase distribution 

A reasonable approximation is that the only contribution to Q* D is from the strength compo- 
nent for D <C D w (the strength regime) and from the gravity component for D 2> D w (the 
gravity regime). Application of the analysis of $23] then suggests that the size distribution 
should have two phases described by power laws (equation [TT} with 

Cfc = (7-a/3)/(2- fl /3), (24) 

in the strength regime and CC;, = (X given by equation fl!4b in the gravity regime. Power law 
indices in this paper, like a a , that refer to equation (11 It are summarised in Table Q] along 
with a summary of diameters that define where these apply. For the dispersal law assumed 
here, a a = 3.63 and a,b = 3.0. 

Although one might expect the distribution to be continuous at D w , closer inspection 
shows that this cannot be the case, since the condition given by equation l|4} should hold 
regardless of the dispersal law, and so should hold on both sides of (and at) the transition. The 
discontinuity at D w was pointed out by O'Brien & Greenberg (2003). The relative heights 
of the two distributions, and so the magnitude of the discontinuity, can be determined by 
including the constants of proportionality in equation (113b - These constants are « K^(l — 

a fl r 1 ( v L/ 2 ) a " /3 2« 1 ~ a ° )/3 for the strength regime, and likewise except with subscript j, for 
the gravity regime, resulting in 

K a _ fT^fvi^-^/ Qt^ 6 ] a ~ 
K b -\il-a b \l) [qI 1 -^ 6 )' } 

Scaling the total mass in the 2-phase distribution to that required, say M tot , means that 
the distribution is now completely defined. In Figure [T] we show with a dashed line the 
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Fig. 1 Steady state size distribution for particles between l^lm and 1000km undergoing catastrophic colli- 
sions. The y-axis is the mass in logarithmic bins centred on sizes given on the x-axis, where bin width is 
given by 5 = 0.01. Three distributions are shown for relative velocities of 0.6, 1.9, 6.5 km s using the same 
law for the dispersal threshold. The distributions are scaled to total masses of 0.01, 0.1 and 1 times the mass 
of the Earth to allow these to be distinguished (though the shape of the distribution is independent of total 
mass). The distributions are calculated using the analytical 2-phase (dashed line) and numerical method 1 
(solid line). 



predicted distributions for the assumed dispersal law for different values of relative velocity. 
Note that we did not change the dispersal law with collision velocity, although this might 
be expected (Benz & Asphaug 1999; Stewart & Leinhardt 2009). The total mass in the 
distributions were scaled to allow them to be distinguished. The main thing to note from 
the analytical distribution (for now) is that the discontinuity is bigger for larger relative 
velocities, as expected from equation (125b . 

It should also be remarked that for low collision velocities, X c can be greater than 1 for 
the largest particles, so that our simplistic view of collisions in £|2.2l implies zero mass loss 
rate for some particles near the top end of the distribution. In practise these particles can still 
lose mass through cratering collisions, illustrating how our assumption that mass loss rates 
scale only with the catastrophic collision rate breaks down in certain regions. Such particles 
could not reach steady state if our assumptions were valid, and a more detailed treatment is 
needed to assess evolution in this regime (e.g., £12. 5b . Rather than assume some primordial 
distribution for such particles, here we set the mass to zero in all bins for which X c > 1 at 
the top end. Note that, although X c can also be greater than 1 at small sizes, this does not 
prevent the destruction of such grains, since they can be destroyed by particles larger than 
themselves. 

2.6.2 Numerical 2-phase with ripples 

There were many approximations in the analytical 2-phase model that are overcome in the 
numerical model. Notably the numerical model should be valid regardless of the shape of 
the dispersal law, and should also be applicable (within the limitations of the assumptions) 
in the situation where the size distribution is finite in extent, being truncated at small and/or 
large sizes. 

Figure[TJalso shows with solid lines the result of numerical method 1 ( £12.4. lb that solves 
equation d 1 5 b - Convergence of this scheme was confirmed by checking that m^ c is indeed 
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the same for all bins, as expected by construction. The first thing to note is the excellent 
agreement with the 2-phase analytical model. There is slight quantitative disagreement on 
the slope of the distributions, which we attribute to the terms that were discarded from 
equation d 1 Ob when deriving equation d!2t . Thus, from Figure[T] we expect the slope of the 
mass distribution in the strength regime to be slightly steeper than predicted by a a , and for 
the distribution in the gravity regime to be slightly shallower than that predicted by (Xb- 

As expected for a realistic distribution, there is no discontinuity at D w . However, as 
discussed by O'Brien & Greenberg (2003, see their Fig. 3) and previously noted by Durda et 
al. (1998), the transition at D w does cause a ripple that extends to larger sizes, and is captured 
by the numerical model. The transition does not affect the distribution in the strength regime 
because, as discussed in O'Brien & Greenberg (2003), the distribution in the strength regime 
is independent of that in the gravity regime, which follows in our analysis because m^ c 
depends almost entirely on the number of smaller objects in the distribution with size around 
X c Dii. We find that the magnitude of the ripple (i.e., the departure from the power law slope) 
depends on the collision velocity, with higher velocities causing larger ripples. The ripple 
decreases in magnitude to larger sizes, with a peak-to-peak wavelength that also depends 
on collision velocity and varies with size (decreasing wavelength to larger sizes). All such 
effects follow naturally from the dependence of X c with size, which is the parameter that 
sets the magnitude and wavelength of the ripple, since for particles of a given size the larger 
X c is the nearer in size particles have to be to affect their number. In the gravity regime X c 
increases with size (and velocity) causing the observed dependence. This also means there 
should be a dependence on Q* D in that weaker planetesimals should also result in a more 
pronounced ripple with a longer wavelength. 

A similar ripple is also seen at small sizes, caused by the fact that the distribution is 
truncated at 1 pm. While this truncation is implemented simply by not considering bins < 1 
pm, a truncation is physically motivated by radiation pressure, since grains smaller than 
Dbi, where 

D w «(2300/p)(L*/M*) (26) 

in pm (for p in kg m~ 3 , and L± and M± in Solar units), are placed on hyperbolic trajec- 
tories as soon as they are created (e.g., Burns et al. 1979). In this case the lack of < 1pm 
impactors causes the equilibrium number of 1pm grains to be enhanced. The consequence 
of this is an enhanced destruction rate of objects of sizes typically destroyed by 1pm grains 
(i.e., those for which X C D = 1pm). This causes a dearth of such grains and so a higher equi- 
librium number of grains that would have been destroyed by them, and so on. Such ripples 
have previously been discussed by Campo-Bagatin et al. (1994) and seen in the simulations 
(Thebault et al. 2003; Krivov et al. 2006). There is a similar dependence on collision veloc- 
ity of the magnitude of the ripple, and its wavelength, in that both increase with v re /. This 
can again be understood in terms of X c , which is both smaller for larger velocities (and so 
particles can be affected by those at much smaller relative sizes), and which decreases with 
size in the strength regime. This ripple would also be expected to be more pronounced for 
weaker planetesimals. The details of this ripple should, however, not be overinterpreted. For 
example, if the truncation is caused by radiation pressure, then grains just above the blow- 
out limit (equation [261 should also have a higher eccentricity than larger particles, which 
numerical simulations show would affect the shape of this ripple (Krivov, Mann & Krivova 
2000; Thebault & Augereau 2007). The ripple may also be less prominent if particles have 
a range of compositions and so a range of truncation sizes and strengths. 

In short, all of the features of the steady state size distribution arising from collisions 
that have been discussed in the literature (power laws and ripples) are reproduced with our 
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simple numerical scheme. This allows the shape of the distribution to be analysed without 
having to evolve size distributions forward in time until they reach steady state, which can 
be time-consuming for large particles with low collision rates, if timesteps are set by the 
small particles that evolve much faster. This should allow us to consider how the steady 
state distribution depends on parameters such as the shape of the dispersal threshold law 
and the collision velocity. But more importantly for this paper, allows us to consider how 
that shape is affected by loss processes. 

2.7 Timescale to damp perturbations from steady state 

The steady state size distributions presented in the preceding sections have some issues at 
the top end of the distribution. It was already discussed that the numerical scheme employed, 
to keep mass loss rate independent of size, should be valid for most of the size range, but 
that for a distribution that is truncated at sizes above D\ this should be invalid for the size 
range in which the redistribution function means that objects of size D\ (and larger, had 
they existed) deposit a significant fraction of their mass. Although a numerical scheme was 
outlined that should be able to account for this (equationll6ll. we should consider the validity 
of the assumption that the largest particles are in steady state. 

First consider a bin to have a mass that departs from that expected when the distribution 
is in quasi-steady state by a small fraction e k , so that 

m k =m ks {\+e k ). (27) 

Since rht° and R c k are independent of for small enough bins (or a small enough perturba- 
tion), substituting equation d27l > into equation ([TJl gives 

e k /e k = -R c k . (28) 

In other words, perturbations are damped on the collision timescale. 

2.8 Analytical 3-phase distribution 

Another way of interpreting equation (I28t is that the collision timescale also defines the 
timescale for a given bin size to reach steady state. Thus for a given evolution timescale, 
t age , all particles with collision rates higher than 1 /t age (i.e., all of the smallest particles) 
would be expected to have reached steady state, which we denote with subscript ss, whereas 
those with smaller collision rates (i.e., the largest particles) would be expected to retain their 
primordial size distribution, which we define here using equation ill 1 t i with a subscript p. 
We define D t to denote those particles with index k = i, for which 

R c k {D,) = \/t age . (29) 

This motivated Lohne et al. (2008) to define a 3-phase size distribution comprised of the 
2 phases in steady state given in £12.6.11 in addition to the primordial phase. These authors 
connected the size distribution in each phase continuously. However, if it is assumed that 
there are two distinct portions of the size distribution, primordial and steady state, with the 
destruction of primordial particles feeding those in steady state (in a manner similar to the 
approach of Dominik & Decin 2003), then there may be a discontinuity at this transition. 
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The reason is that bins in the steady state population have mass input equal to that 
lost. This means that the rate of mass loss from the steady state population is the sum over 
this population considering just the mass put into particles small enough to fall outside the 
cascade 

hum 00 

™.7 s = c£ £ F(k-i), (30) 

i=i, k=i mi „ 

where we have used the fact that m ( ^ e = C in the regime considered to take it outside of the 
sum, and i mm = for initial considerations in this paper, but would be the larger of z'^ and 
i pr if a consideration of P-R drag was necessary ([J3jl. Since C scales with the square of the 
total mass in steady state (for fixed i,), so does the mass loss rate. The mass input to the 
collisional cascade, however, is set by the sum over all bins in the primordial distribution, 
considering just the mass put into sizes in the collisional cascade 

H hiiin 

!=1 k=i, 

where rhj c is set mostly by the primordial distribution. 

Thus in this approximation, the mass in the steady state population represents a balance 
between the two, i.e., m+ = rhj s so that 

C = k ~" (32) 

m,L^i mia F(k-i) 

and does not necessarily lead to a continuous distribution at D,. This is particularly so if D, is 
close to D\, since at this point the mass input to the steady state population from primordial 
particles is low, and so that population must also have low mass and the discontinuity is 
large. 

In fact the transition from primordial to steady state is not expected to be discontinuous, 
and such a discontinuity may compromise the validity of the scale independence assumption 
of our redistribution function ( i]2.5b . Nevertheless it is hoped that this approximation still 
captures the important physics. For example, although the discontinuity is large in the regime 
in which D, approaches D\ , the primordial population may, in some respects such as the long 
term evolution of the steady state population, be irrelevant at this point. This is because, 
in £12.91 the system age is set to be equal to the collisional lifetime of objects of size D t , 
which may well have a negligible contribution from the primordial population. If this is the 
case, then the collisional lifetime would be expected to scale inversely with the mass in the 
steady state population, leading to disk mass falling off inversely with age, as expected from 
previous analytical and numerical studies (e.g., Dominik & Decin 2003; Wyatt et al. 2007a; 
Lohne et al. 2008). 



2.9 Evolving the size distribution 

The scheme of £12.81 allows us, for a given primordial distribution and transition size, to 
define the discontinuity and so derive a 3 phase distribution analytically. Further equating the 
collision lifetime of objects of size D, with the system age would also permit the evolution 
of an input primordial size distribution to be determined, in a similar manner to Lohne et 
al. (2008). Here we follow a comparable approach, except that the 3-phase distribution is 
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derived using numerical method 2 ( £j2.4.2t i.e., using equation IT6ll. with the largest size 
bins (above D t ) fixed at the primordial size distribution. This automatically captures the 
discontinuity (and its effect on the collision rate of particles at the transition) and achieves 
steady state in the appropriate size range. 

To consider how accurately our scheme reflects simulations that solve for the evolution 
using the full equation (fT) (or its equivalent), we compare our results with those of the 
ACE ran ii-0.3 of Lohne et al. (2008). That run considered the evolution of a 7.5-15AU 
planetesimal belt with a maximum eccentricity of 0.3, maximum inclination of 0.15 rad, 
and largest planetesimals 148km in diameter, for a Sun-like star. The primordial distribution 
had 1M® distributed according to equation i ll 1 1 with a p = 3.61. The result of our scheme 
for these parameters, with 8 = 0.02 (and so 1279 diameter bins), is shown in Figure [2] 
Figure [2] should only be compared qualitatively with figures in Lohne et al. (2008), since 
that paper mistakenly used a slightly different gfo = 0.014 J kg~' , a factor of 2.5 higher than 
stated (and used in our nominal runs, £|2.6b (T. Lohne, priv. comm.). However, a quantitative 
comparison is given within Figures |2j;,|2]i and[5£, that also show ACE results for run ii-0.3 
with Qj, = 5.6 x 10~ 3 J kg -1 (A. Krivov, priv. comm.). This comparison shows that our 
scheme mostly reproduces the same features, both qualitatively and quantitatively, but there 
are also some qualitative and quantitative differences that merit further discussion. 

For example, the size distribution at specific ages (Figure |2^) compares well with that 
of Lohne et al. (2008). Although the details of the ripples close to the blow-out limit are 
not identical, they are qualitatively similar, and differences are expected due to the models' 
different assumptions about collisional outcome and in their treatment of the dynamics of 
grains close to this limit. However, the discontinuity at the transition size is not apparent in 
the Lohne et al. simulations. Figure^ illustrates why this discontinuity cannot be physical, 
since although the mass evolution of different bins shows generally good agreement, all 
sizes undergo a sudden loss of mass at the time they are assumed to reach steady state, 
whereas mass evolution should be continuous. For this reason we do not expect our scheme 
to reproduce details of the size distribution at the transition between primordial and steady 
state, such as the divot found by Fraser (2009). 

The most instructive plots for comparative purposes are those for the evolution of total 
mass (Figure |2];) and dust mass (Figure [2jl). Again, these show qualitatively good agree- 
ment. For example, total mass exhibits a slow fall-off with time from O.lMyr, but speeds up 
around lOOMyr tending to a fall-off with age °c r~ 0,94 , close to the inverse fall-off with time 
predicted above. However, the total mass lost at 20Gyr is ~ 92% of the initial mass, whereas 
the ACE ran loses closer to 81% on the same timescale, indicating that mass loss is quan- 
titatively faster in our scheme. The dust mass shows similarly good qualitative agreement, 
with dust mass falling off between l-200Myr °c t ~ 0,34 for our scheme and °c t ~ ' 41 for ACE. 
Although it is noticeable that our scheme underpredicts the dust mass by a factor ~ 2 at 
most ages, this is a relatively small difference given the many orders of magnitude in mass 
over which the size distribution is being considered. However, our scheme also predicts a 
turnover beyond ~ 2 Gyr, tending to dust mass falling off °c r - 91 . This is in qualitative 
agreement with Lohne et al. who predict an inverse fall-off with age at late times, once the 
largest planetesimals have come to collisional equilibrium, except that the turnover in their 
plot is not noticeable since it occurs much later. 

Figure [2£ shows that the dependence of dust mass on total mass in our scheme is very 
similar to that of Lohne et al. (2008), apart from the small underprediction of dust mass 
noted above. This suggests that a reinterpretation of age could lead to greater agreement 
between the schemes. Remember that the size distribution is worked out for a given D t , 
and the corresponding age is then derived by working out the collision rate of objects at 
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Fig. 2 Evolution of a planetesimal belt with the parameters quoted for ACE ran ii-0.3 of Lohne et al. (2008). 
For comparison, the results for such an ACE ran are reproduced on (c), (d) and (f) with a dotted line (A. 
Krivov, priv. comm.), whereas for (a) and (b) qualitative comparison with figures in Lohne et al. (2008) is 
possible (see text for details), (a) Mass (m^) in different size bins (D^) at times of 0.05, 0.5, 5, 50, 500 Myr and 
5, 20Gyr (compare with Lohne's Fig. 4 bottom), (b) Mass (m^) in 50 different (but equally-logarithmically- 
spaced) size bins as a function of time (compare with Lohne's Fig. 4 top), where bin size could be deduced 
from plot (a), (c) Total mass as a function of age. (d) Dust mass (in particles D < 2mm) as a function of 
age. (e) Dust mass as a function of total mass, (f) Transition diameter (from primordial to steady state) as a 
function of age. In practise, the size distribution is calculated for a given D, , which is then converted into an 
age by determining the collision lifetime of objects at the transition. 
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the transition and equating that with the inverse of the age (see Figure[2j). In Figure [2]we 
consider some of the model parameters that may affect evolutionary timescales. First, since 
collisions with the primordial population may shorten the collision lifetime at late times 
(when such a population is not physically realistic), the dashed line shows what happens 
when collision rates are calculated assuming that steady state and primordial particles can 
only collide with other members of their respective sub-population. This is implemented by 
changing the indices over which the collision rate is summed in equation Q. As expected 
this results in slightly longer evolutionary timescales for both total and dust mass. Notably 
the evolution is now closer to that found by the ACE simulations in that 85% of the total 
mass has been lost at 20Gyr, and the turnover in the dust mass evolution at 2Gyr is no longer 
present. We also considered the effect of the redistribution function, by trying a cratering- 
like function weighted toward objects of size (with r\ rmax = 1 and a r = 3.9, see dash- 
triple-dot line). This has very little effect on the evolution of either total mass or dust mass, 
and it can be concluded that the size distribution in the collisional regime has very little 
dependence on the redistribution function (but can be used as a probe of the latter in the P-R 
drag regime, ij3). 

Thus, despite this limited comparison, it seems that our prescription matches the results 
of more detailed simulations without the need for such time-consuming simulations to deter- 
mine the evolution on long (>Gyr) timescales. Although the 3-phase analytical prescription 
given in Lohne et al. (2008) does likewise, our model has the benefit that we reproduce 
more detailed features of the size distribution in the steady state regime (such as ripples), 
and moreover it is possible to consider additional effects such as loss processes like P-R 
drag (see Sj3j. We expect that our scheme would replicate the results of other sets of full 
simulations (e.g., with different assumptions about collision outcome), though a detailed 
comparison would be needed to quantify (and improve through refining the implementation 
of the scheme) the level of agreement. 



3 Steady state with P-R drag 

Despite the uncertainties in the details of how the size distribution transitions from primor- 
dial to steady state, equation d28b demonstrates that the small size end of the distribution is 
firmly in steady state and so can be treated as such. Consider now a belt in which both col- 
lisions and loss processes are acting. In this section we will consider the effect of Poynting- 
Robertson (P-R) drag. Poynting-Robertson drag is the dominant loss process for dust in the 
asteroid belt. 



3.1 Analytical fourth phase 

To work out the loss rate from the belt we consider the timescale for dust particles to migrate 
from r ml d to r m!£ f — dr/2. Clearly, such particles have not yet been removed from the system, 
but are then fed into an inner region which should then be considered in some other way; 
our prescription means that we are only considering the size distribution of particles in the 
belt region. The loss rate by P-R drag from the belt for individual particles is 



R{ — A pr D k 1 , 

A pr = 2.9 x \{r 6 L k l{M*Pdrr^ d {\-0.7Sdr/r nM )), 



(33) 
(34) 
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where r ml d and dr are in AU, L* and M* in solar units, and p in kg m 3 , to get a collision 
rate in yr _1 with D k in m. Thus the loss rate from the k-th bin is 

m- pr = m k A pr D- 1 . (35) 

Because m k « DZ~ a (from equationsl9land[TT1>. then in the steady state coUisional regime 
where m k R c k is independent of D k , the loss rate from collisions is 

Thus in both the strength and gravity regimes this drops slower than P-R drag, which falls off 
as R p k ' cx D k l . Thus it is possible that at some size, which we call D pr (with a corresponding 
index i pr ), these two are equal, and so 

R c k (D pr )=R pr {D pr ). (37) 

For sizes smaller than D pr the loss is dominated by P-R drag, whereas for sizes larger than 
D pr losses are dominated by collisions, an assumption that will be justified later on. For some 
disks D pr is smaller than the size for which particles are removed by radiation pressure (Dm), 
and the mass of particles removed by P-R drag is insignificant. In fact this is necessarily the 
case in most debris disks that are bright enough to be detectable (Wyatt 2005). However, 
it is not the case for the Solar System's zodiacal cloud, and it will not be the case for the 
expected population of fainter extrasolar debris disks we have yet to discover. 

First let us work out what the steady state distribution is for the size range in which P-R 
drag is the dominant removal mechanism. If we assume that the rate of coUisional loss in this 
regime m7 c <C m7 pr and so can be ignored, then in steady state the equivalent to equation 
lO is 

^ pr ='tmr/F(k-i), (38) 

where the sum only extends up to i = i pr , since mass loss from i > i pr is assumed not to be 
coUisional. Since m7 c = Q is approximately independent of size in the coUisional regime 
(i.e., C k — C), this can be taken outside the sum giving the steady state distribution as 

m ks =A pr l CD k ^F(k-i). (39) 
i=l 

Thus the size distribution in the P-R drag regime has a slope that is set by the redistribu- 
tion function modified by an additional factor of D k . This prescription for the distribution is 
continuous at D pr because at that location the sum in equation l !39l > is unity and so the mass 
calculated in the P-R drag regime at that location is m ks (D pr ) = C/R k " , which is the same 
as that calculated in the coUisional regime m ks (D pl .) = C/R c k . For the redistribution function 

defined in £|2.5l bv a r , the sum works out to be F{k - i) cx (D pr /D k ) a '- 4 , which means 
that 

m ks *A- r l CD pr {D k /D pr ) 5 - a '. (40) 
Or in terms of equation (II It with subscript pr, 

a pr = a r -l. (41) 

This is important because it means that the size distribution of small particles in the P-R 
drag regime is an indirect measure of the slope in the redistribution function. 
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Since a r < 4 the size distribution in the P-R drag regime is depleted in small particles. 
The main consequence of this is that the rate of collisional losses decreases toward smaller 
particles. This is because the mass loss rate due to collisions given by equation d 13b applies 
regardless of whether the slope in the size distribution a is set to make the exponent in 
this equation equal to zero (as is the case in the collisional steady state regime). Given that 
m k °c D A k ~ a this means that, for particles that are small enough to be in the strength regime, 

4ocD^"- (1 - a)fl/3 . (42) 
So in the P-R drag regime where a is set by the redistribution function, for the parameters 

3 8 9ff 

assumed here, R c k °c Tn us, as we have assumed that a r < 4, collisional loss rates 

increase with size in the P-R drag regime, meaning that collisional losses can be ignored for 
all particles smaller than D pr , justifying the earlier assumption. 

3.2 Numerical size distribution 

The steady state distribution, including the influence of P-R drag, can be found numerically 
by extension of method 2 ([J2A2}. Again we equate mass loss with mass gain in each bin, 
to derive a revised version of equation H6) 

m ks =C k /{R c k +A pr D- y ), (43) 

where Q is again given by equation dl7l i. 

This method is used in Figure [3] to illustrate the points raised in i)3.1l regarding the 
resulting distribution. To simplify the discussion, these simulations use exactly the same 
parameters as those used in the comparison with run ii-0.3 of Lbhne et al. (2008) that was 
discussed in £|2.9l and Figure[2] The main difference is that here we include the effect of P-R 
drag, and we also vary the parameter cc r . Since the effect of P-R drag only becomes notice- 
able with these parameters for ages ^> lOGyr, its exclusion from the previous section (and 
from Lohne et al.'s simulations) does not compromise those results. However, this means 
that we are considering here the distribution at ~ 500 Gyr ages that are greater than the 
age of the Universe. Although it seems ridiculous to consider such ages, this is not a major 
concern, since the same features would appear on shorter timescales in belts with different 
parameters (e.g., with smaller D\ or r m! y), or in those started with low total mass, and effects 
such as dynamical depletion can be modelled using collisional evolution timescales longer 
than the system age (e.g., Bottke et al. 2005; 

Figure [3^ shows how P-R drag causes a turnover in the size distribution toward smaller 
particles, and how the slope in the P-R drag regime is set by the slope of the redistribution 
function, which is exactly as expected (remembering that the slope on this figure is expected 
to be 5 — cc r ). It also shows how insensitive the rest of the size distribution is to the redis- 
tribution function. It is evident that the ripples in the size distribution at small sizes have 
disappeared, which is expected above the blow-out limit, as the loss of particles just above 
this limit is no longer controlled by collisions. However, the turnover at D pr has the potential 
to cause a ripple in the size distribution just above this size, which is not seen presumably 
because the turnover is a gradual rather than sharp feature. 

As for the location of the turnover, Figure |3j> shows that the diameter for which R c k = R? r 
is D pr = 660, 430, 260 /im for a r = 3.0, 3.5, 3.9, respectively, noting that R c k is calculated 
from the size distribution of Figure [3^ rather than that with no consideration of P-R drag 
(which from extrapolation of the power law in the collisional regime would have resulted in 
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Fig. 3 Steady state size distribution of a planetesimal belt with the same parameters as assumed in run 
ii-0.3 of Lohne et al. (2008), but including loss due to P-R drag. Different assumptions about the slope 
in the redistribution function, a r = 3.0,3.5,3.9 are shown on all plots with dotted, solid and dashed lines, 
respectively. The evolutionary age was set to give the same mass of lm particles for the different assumptions, 
corresponding to an evolutionary age of 290, 540 and 2100 Gyr for the different a r respectively, (a) Mass 
(m k ) in different size bins (D k ). (b) Loss timescales for particles of different sizes in these distributions due 
to collisions and P-R drag, (c) Mass loss rate from the inner edge of the belt due to P-R drag. The dash- 
dot line shows the polynomial fit to the accretion rate of dust onto the Earth from the LDEF satellite (Love 
& Brownlee 1993) scaled up by a factor 1 / (40D k ) to account approximately for the accretion efficiency to 
recover the size distribution in the vicinity of the Earth. 
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a smaller value of D pr sa 150 Jim that is independent of a r ). For all values of a r , this is a 
factor of 1.9-2.0 times larger than the sizes obtained from the intersections of power law size 
distributions describing the P-R drag and collisional regimes, and a factor of 1.5-10 times 
larger than the peaks in the cross-sectional area distributions. 

The final Figure |3j: shows the mass loss from the belt due to P-R drag. Since this mass 
can be considered as a gain term for a population of particles that are evolving due to P-R 
drag we denote this as m£ pr which is °= m.]fl~^ , and so from equations ((9} and (lilt is also 
°c Dl~ a . Thus in the P-R drag regime, D < D pr , this distribution is set by the redistribution 
function 

™ + k pr {D<D P r)~D A - a ', (44) 
whereas in the collisional regime, D > D pr , it is set by the strength law 

mt pr (D>D pr )o,D ( f 3 - 2a ^ 6 - a \ (45) 

As mentioned above, this mass is not lost from the system, but continues to evolve toward 
the star due to P-R drag, also suffering collisions on the way. It is beyond the scope of 
this paper to include the extra radial dimension needed to consider the evolution of this 
material properly. However, we can consider what happens if collisions no longer amend 
the size distribution, and the particles also pass a planet onto which they may be accreted. 
In this case the size distribution of accreted particles would simply be given by Figure [5J; 
modified by a size dependent accretion efficiency P acc . Accretion efficiencies may plausibly 
vary P acc oc D k due to the slower motion of larger particles past the planet (e.g., equation B2 
of Wyatt et al. 1999), though the detailed dynamics may be more complicated (Kortenkamp 
& Dermott 1998). 

The assumption that collisions can be ignored is likely too crude to be useful for specific 
applications; in particular, collision timescales may be expected to be shorter than P-R drag 
timescales for particles larger than D pr . However, it is notable that this size distribution is 
not so dissimilar to that derived for the accretion of dust onto the Earth by LDEF (Love & 
Brownlee 1993), a scaled version of which is also shown on Figure [3};. This suggests that, 
with further work, it may be possible to leam about the redistribution function from such ob- 
servations. For example, if the collisionless and P acc <x assumptions held, then Figure[3j: 
shows that the LDEF accretion rate for < 200 ^xm particles would imply a redistribution 
function with a r slightly below 3.5, perhaps with a turnover to a smaller value below a few 
10s of jiim. Griin et al. (1985) also noted that the size distribution in the zodiacal cloud at 
1AU may be indicative of the redistribution function, from a consideration of collisional 
gains and losses due to collisions and P-R drag, though these authors favoured the inter- 
pretation that the zodiacal cloud is not in steady state to altering the redistribution function 
from their assumed value of a r = 3.5. 

In fact the redistribution function is poorly constrained. This distribution for the break- 
up of large objects can be measured from asteroid families and is known to depend on impact 
energy, partly due to geometrical considerations (Tanga et al. 1999), but is typically found to 
range from a r = 3.5 — 4.5. Collisional experiments of (smaller) rocks (Fujiwara et al. 1989) 
give distributions that have a r = 2.8 — 5.8 (though 3.7-4.0 is more common). However, it 
is noted that the distributions can turn over to lower values of a r < 3.5 at fragment sizes 
smaller than 1mm (e.g., Flynn et al. 2009), i.e., in exactly the size range of interest, also 
changing to even lower values for < 100 pLm. Thus, although the distribution implied by 
Figure [3}; would appear to have a relatively small value of a r , this is not without precedent 
within the context of previous studies, especially for this size range. 



20 



M.C.Wyatt et al. 



3.3 Evolution in P-R regime 

The evolution in the P-R drag regime was discussed in Dominik & Decin (2003) who found 
that the dust luminosity should decrease <x \/t 2 in this regime. We find that the situation 
is substantially more complicated than the prescription given by these authors, since we 
treat an evolving size distribution in which the numbers of particles of different sizes evolve 
differently. Nevertheless, for certain assumptions we recover the same dependence for the 
evolution of luminosity. 

To consider the evolution in the P-R drag regime, we first need to determine the evolution 
of D pr . Here we will consider the situation in which the largest particles are in equilibrium, 
so that the mass in all bins above D pr are being depleted with age as l/t. The implications for 
Figure(3j) are that the collision lifetimes of all particles (above D pr ) increase <x t, and so the 
turnover evolves to larger sizes. Here we ignore the effect of the turnover on the calculation 
of the collision time at D pr ; this means that we miscalculate D pr , but by a fraction that 
should be the same at all ages for a given a r , and so should not affect conclusions about the 
evolution. Thus equating collision and P-R loss rates we find that 

D (lr cxf(M/P+2«) i (46) 

or D pr °c t 158 for the parameters of our simulation. 

For particles larger than D pr , their number continues to deplete inversely with age. How- 
ever, for particles that are now in the P-R drag regime, the increase of D pr with age means 
that they are being replenished by a population for which the size at which its bottom end is 
truncated is also increasing with age. The number of particles in the P-R drag regime thus 
depletes much faster 

m ks {D <D pr )oc t [«r(6-«)-30]/(3+2«) ; (4?) 

or °c r~ 2 79 for the parameters of our simulation with a r = 3.5. The evolution that is observed 
then depends on what particle sizes are being probed. For example, the evolution of dust 
mass defined as the mass in D < 2mm particles would be expected to fall off inversely with 
age at early times when D pr <C 2mm, but tend to a fall off °c ;~ 2 - 79 by the time D pr 3> 2mm, 
and likewise for similar definitions of dust mass. 

However, for certain observations the contribution to the emission from the disk is dom- 
inated by the total cross-sectional area in the distribution. For the distributions considered 
here, with a r < 4, that area is dominated by particles of size ~ D pr . Although this size 
evolves in a way that depends on the strength law (equationl46ll. it turns out that the evolution 
of total cross-sectional area does not, since a lol m^(D re f)D p 3 ,^ 2a ^ i ' a 6 \ where m^(D re f) 
is the mass in particles at some fixed reference diameter D re f that is above D pr . In the 
regime we are considering in which the largest particles have reached collisional equilib- 
rium m\SPref) x t _1 , and so from equation ( 146 1 we find that o tot <x t~ 2 as predicted by 
Dominik & Decin (2003). Neverthless, all observations are prone to bias toward probing 
certain particle sizes (see, e.g., Fig. 5 of Wyatt & Dent 2002), and so one would expect all 
observations to tend to fall off at a rate given by equation (1471 ) on long enough timescales, 
and it is important to consider the particle sizes contributing to a particular observation to 
understand how disk brightness is expected to evolve. 

The above points are illustrated in Figure |4] which shows the evolution from an age of 
lOGyr of the planetesimal belt considered in Figure [2] but including P-R drag at all ages 
(note that a r = 3.5 in these simulations). Figure^ shows how the size distribution at lOGyr 
is unaffected by P-R drag (justifying its exclusion from the simulations shown in Figure [2}, 
but that soon after this the turnover due to P-R drag appears at a size that gets larger with 
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Fig. 4 Evolution of the planetesimal belt considered in Figure [2] on timescales long enough for P-R drag 
to affect the evolution of the dust. Although these timescales are unrealistically long for this parameter set, 
shorter timescales may apply for different belt parameters, and dynamical depletion can result in evolutionary 
timescales that are longer than the system age ( £|4. 1 .2i . (a) Mass (m^) in different size bins (Z)j) at ages of 10, 
30, 100, 300, 1000, 3000, 10000 Gyr. (b) Loss timescales for particles of different sizes at these ages due to 
collisions and P-R drag, (c) Evolution of total mass, dust mass (both defined as mass in particles smaller than 
2 mm, and smaller than 10 lim), and cross-sectional area. 
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Table 2 Evolutionary stages for the belt discussed in Figs.[2]and[4] with the ages over which the stage lasted, 
the phases present in the size distribution (P=Primordial, S=Strength regime, G=Gravity regime, PR=P-R 
drag regime), and the dependences on time of total disc mass, dust mass and total cross-sectional area. 



Stage 


Ages 


Phases 


m ta t 


™dust 










P 








I 


< 0.5Myr 


S-P 


flat 


flat 


flat 


II 


0.5Myr-2Gyr 


S-G-P 


gradual 


oc f-0-34 


oc f-OM 


III 


2 - 30Gyr 


S-G 


«f-l 


«r' 


=cr-' 


IV 


> 30Gyr 


PR-S-G 


ocr- 1 


oc f-2-8 


oc ( -2 



age. It is also evident that the ripple above D pr that was noted to be absent from Figure [3]is 
now apparent at late ages. This is because at larger values of D pr , particles at this transition 
have lower values of X c , i.e., they are destroyed by particles that are much smaller relative 
to themselves than those with smaller values of D pr , which means that the turnover appears 
more like a truncation to them. Figure |4j> shows how the turnover size is indeed set at 
all ages by the location where collision and P-R drag timescales are equal. This justifies 
the analytical work in the last few paragraphs, which also provides excellent quantitative 
agreement with the evolution of mass and cross-sectional area seen in Figure [4}:; indeed at 
late times the simulations find m tot <x r~ 10 , rtidm x t~ 2,& , O wl °= r~ 2 - . 

A word of caution is necessary on the evolution implied by equations d46b and J47K 
which is that these apply when the largest objects in the distribution are in collisional equi- 
librium and so are being depleted <x 1 jt. This may not be the case for belts that were simply 
born with low mass, or have undergone dynamical depletion (like the asteroid belt), or have 
enhanced levels of drag (e.g. £|4.3l l. For such systems it would be necessary to consider the 
rate of depletion of particles with sizes above D pr , and substitute that into the above analysis 
in place of the <x \/t dependence assumed there. 



3.4 Summary of evolutionary stages 

To summarise, if we ignore the ripples in the distribution, then the evolution of the size 
distribution shown by the runs of Figs.|2]and|4]show four main stages of evolution that are 
summarised in Table [2] 

The evolution of planetesimal belts with different initial parameters would be expected 
to show similar behaviour. However, there would be quantitative differences. For example, 
the ages of the transitions between the stages would depend on the initial parameters, since 
they correspond (e.g., Table [T} to the collisional lifetime of objects of size D w (transition 
I-II), and of size D\ (transition ITIII), and the age at which the collisional lifetime of objects 
of size Dm is equal to their P-R drag lifetime (transition III-IV). The fall-off of dust mass 
and area with age during stage II also depends on the primordial size distribution and on 
the dispersal law in the gravity regime (see equation 43 of Lohne et al. 2008), and the dust 
mass in stage IV depends on the dispersal law in the strength regime and the redistribution 
function (equation I47t. There may also be qualitative differences. For example, if the pri- 
mordial distribution contained no particles in the gravity regime, i.e., if D\ < D w , then there 
could be no "G" phase in the size distribution, and stage II would be skipped. Similarly, if 
the primordial distribution was made up of different power laws, then there would be differ- 
ent phases associated with the timescale for collisional equilibrium to reach the transitions 
between the primordial power laws. And it was already noted that it may be possible for P-R 
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drag to become important during stage II in which case an evolutionary stage involving the 
4 phases "PR-S-G-P" would be reached. 

4 Steady state with generalised loss processes 

The method of Sj3]for considering the effect of P-R drag on the steady state size distribution 
of planetesimal belts undergoing a collisional cascade can be applied to other loss processes 
in a similar manner. In particular, if the loss process is described by 

ml 1 = m k R' k , (48) 

then its effect on the steady state size distribution can be determined numerically by extend- 
ing method 2 ( £|2.4.2I > to solve 

m ks = C k /(R c k +R l k ). (49) 

Here we consider several loss processes, and their qualitative effect on the size distribution 
(and its evolution), but defer a detailed treatment to later work. Treatment of other loss 
processes such as Yarkovsky forces and gas drag are also possible within this framework, 
as are erosive loss processes, such as sublimation and sputtering, by invoking the gain term 
m^ g in equation <[TJ) (since loss from one size bin results in gain into another). However, 
such considerations are again beyond the scope of this paper. 

4. 1 Power law loss rate 

Consider the action of a loss process that can be defined by equation ( 148b with 

fl[°cZ\ 7 ', (50) 
which is a generalisation of loss by P-R drag for which 7/ = 1. 

4.1.1 Large 7 

As long as y is large enough, the analytical treatment is identical to that for P-R drag, in 
that there is a transition at the size D\ for which 

4(A) =4(A). (51) 

For D < D/ the removal of particles is dominated by the loss process, while for D > D/ 
removal is dominated by collisions. A consideration of the analogous Figure |3p shows that 
large enough here means 7/ > 7/,,„ where jiim ls determined from a comparison of the slope 
in the loss regime (equation I50t with that in the collisional regime (equation [36}. Assuming 
the transition occurs in the strength regime (i.e., D/ < D w ), equation (124b means that 

7/»»=4-a fl = (3-3a)/(6-fl), (52) 

which for a = 0.3 means 7/ > 0.37. Thus the evolution is qualitatively similar, and the size 
distribution in the loss regime is still indicative of the redistribution function, except that the 
index is now 

a, = a r -y,. (53) 
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The mass loss from the belt by loss processes is riitj « D% m for D > D/, but is still 
rhzf oc D^ a '' for D < D). Assuming the largest particles are in collisional equilibrium so 
that particles in the collisional regime are being depleted <x t~ x , then by equating the rate of 
collisional loss R c k °= Z)™° -4 ; -1 with the rate from loss processes (equation 150). we find that 
the transition size follows 

1 

Dioctn-rim. (54) 

For similar reasoning to that used to derive equation (147 K this means that the mass in small 
particles below the transition evolves 

a r -4-2{y l -r lim ) 

m ks (D < Dfi « t ■n-rum . (55) 

4.7.2 Dynamical depletion 

However, if Yi < Yum, then it is not small particles that are most affected by loss processes, 
but particles above D\. This is the situation for most dynamical loss processes for which 
loss rates are independent of size and so yi = 0. The effect on the size distribution is more 
complicated in this case. 

For example, consider the belt of Figure [2] that has evolved to t c i ep = 0.5Gyr without 
dynamical mass loss, but then undergoes impulsive and significant dynamical depletion by 
a factor fd ep . This would cause the size distribution to be retained, but to drop at all sizes 
by the same fraction. This would mean that additional evolutionary time would be needed 
for the smallest objects in the primordial distribution (~ 27km in diameter) to come to 
collisional equilibrium and continue to evolve. In the meantime the size distribution in the 
portion that was originally in steady state (D < 27km) would be retained since the depletion 
would not have affected the balance of mass gain or loss from such bins. 

This example illustrates two ways in which dynamical depletion can affect the size dis- 
tribution of what appears to be primordial at late times: some parts of the distribution may 
represent a scaled down version of the original primordial distribution, whereas others may 
represent the steady state distribution attained during an earlier high mass phase of evolu- 
tion, even though their collisional lifetime is longer than the current system age. This point 
was recognised by Bottke et al. (2005), who showed how dynamical depletion can be mod- 
elled by considering collisional evolution that takes place over timescales longer than the 
system age. This is in agreement with our analysis, since we expect the size distribution of 
the above example following the depletion to look like that expected for a primordial dis- 
tribution that is scaled down by fd ep then evolved for a time 1 / fdep times longer than the 
current system age. As noted by Bottke et al., it is not possible to infer fd ep from the current 
distribution (of the asteroid belt in their application), i.e., we cannot determine if the original 
distribution was very massive and short-lived, or less massive and longer-lived. However, it 
is possible to model the evolution of the size distribution for depletion by a known amount 
fdep at a known time t& ep . 

4.2 Realistic radiation forces 

The prescription for the effect radiation forces in the preceding sections is an approximation 
of the true physics. For example, the truncation of the size distribution by radiation pressure, 
which occurs for particles with D < is implemented by simply not including such bins 
in the analysis. However, a more realistic approach would be to include radiation pressure as 
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a loss process with an associated (short, i.e., less than orbital) timescale for bins below Dm- 
Thus it would be possible to account for the possibility that the blow-out population collides 
with bound grains (Krivov et al. 2000), and for a gradual rather than sharp truncation due to 
a finite eccentricity of the parent particles (which means that jS > 0.5 particles can remain 
bound) and/or due to a range of particle properties. It was already noted that the increase 
in eccentricity for particles close to Dbi is not accounted for in this model, but again this 
could be incorporated using an appropriately modified intrinsic collision probability in 
equation Q (e.g., Krivov et al. 2005; Wyatt et al. 2010). 

Another approximation is the inverse linear dependence of the P-R drag loss rate on par- 
ticle size. In fact R? r °c /3, and although /3 °= l/D for large black body particles, this function 
peaks at 0.1-1 jiim particles (i.e., at particle sizes that are comparable to the dominant wave- 
length of stellar radiation), turning over for smaller particles and tending to jS independent 
of size for particles <0.1 /im (e.g., Gustafson 1994). For particles around high luminosity 
stars this turnover is inconsequential, since it occurs at sizes below D/,/. However, for low 
luminosity stars it is so important that there may be no particles that can be removed by 
radiation pressure, i.e., /3 < 0.5 for all particle sizes (e.g., Sheret et al. 2004). The conse- 
quence of this is that, in the absence of other loss processes, P-R drag is the only removal 
mechanism for dust from the cascade, and so it is important to consider the effect of the 
turnover in j6 vs D on the steady state size distribution. A reconsideration of equation d40t 
shows that the size distribution in this regime should really be 

m ks ~D 4 - a '/l5. (56) 

This goes some way to explaining the shape of the size distribution seen in Fig. 4 of Rei- 
demeister et al. (2010), since equation d56t would predict a distribution of drag dominated 
particles that is fairly flat in the 0.04-10 |im size range with a dip at ~ 0.4 jiim. 

4.3 Stellar wind 

Although the discussion has focussed on radiation forces, it is notable that P-R drag only be- 
comes important for debris disks that are so low in density that they are no longer detectable 
with current instrumentation (Dominik & Decin 2003; Wyatt 2005). However, stellar wind 
forces act in exactly the same way as radiation forces. Their pressure component (which is 
analogous to radiation pressure) can usually be ignored, but their drag component (which 
is analogous to P-R drag) can significantly exceed P-R forces (e.g., Plavchan et al. 2005; 
Reidemeister et al. 2010). Thus stellar wind drag can be considered as an enhanced form 
of P-R drag using the same analysis but with a constant A sw in equation d33 b that means its 
effect on the size distribution can become apparent at much younger evolutionary ages (i.e., 
for much higher, and potentially detectable, disc masses) than its radiation counterpart. As 
per the discussion of £14.21 it may also be the dominant loss mechanism for the belts of low 
luminosity stars. No change in the qualitative picture of disc evolution presented in i)3.4l is 
necessary, except to note that the turnover in the size distribution may become apparent in 
earlier evolutionary phases (e.g., in phase II as already noted in ^3 -4b - 

5 Conclusion 

In this paper we presented a new scheme for determining the shape of the size distribution, 
and its evolution, for collisional cascades of planetesimals undergoing destructive collisions 



26 



M.C.Wyatt et al. 



and loss processes like Poynting-Robertson drag. This scheme, outlined in 32 considers 
mass gain and loss in different size bins, and treats the steady state portion of the cascade by 
solving for the assumption that = 0, and assumes that larger objects retain their primor- 
dial distribution. Many results can be obtained analytically, providing insight into the origin 
of the shape of the size distribution. For example, for particles in the steady state collisional 
regime this prescription was shown to lead to mass loss rate being independent of size in log- 
arithmic size bins; thus the steady state size distribution is independent of mass gain into the 
bins in this regime. However, perhaps most notably this scheme provides an efficient and 
readily implemented numerical method for calculating the size distribution. Although the 
scheme as presented assumes a simple scale independent redistribution function that defines 
the average mass distribution of fragments produced in collisional destruction of objects of 
a given size, and considers only the size distribution in one radial bin, it may be expanded 
to more complex redistribution functions, and to include a radial dimension (or to include 
dimensions of orbital elements). 

Our scheme, and its application to consider the evolution of the size distribution in the 
collision-dominated regime, was shown to reproduce qualitatively all of the features found 
by more detailed (and computationally more time-consuming) models that step size distri- 
butions forward in time. Such features in the size distribution include ripples above the ra- 
diation pressure blow-out limit, the change in slope at the transition from strength to gravity 
regime, the ripple above this transition, and the transition to primordial regime. For evo- 
lution, the scheme reproduces the initially slow fall-off of dust mass with time °c f~ - 34 
followed by steeper decline °c t , as well as an evolution of total disk mass that is much 
flatter in the initial phases but also tends to °= t _1 at late ages. 

One of the main purposes of the paper was to consider the effect of loss processes on 
the steady state size distribution, which was done in detail for P-R drag in |J3] This was 
found, both analytically and confirmed numerically, to cause a turnover at small sizes to a 
size distribution that is set by the redistribution function. This means that information about 
the redistribution function may be recovered by measuring the size distribution of particles 
undergoing loss by P-R drag. This was illustrated by comparison with the size distribution 
of particles accreted by the Earth, though it was noted that further work on the collisional 
evolution of material inside the source belt and on the accretion efficiency as a function of 
particle size would be needed before definitive conclusions on the redistribution function 
could be reached. The evolution of the size distribution when P-R drag is included is sum- 
marised in £]3.4l and Table[2] Notably it is found that although cross-sectional area drops with 
age cx f~ 2 in the PR-dominated regime, dust mass falls <x f~ 2 - 8 , underlining the importance 
of understanding which particle sizes contribute to an observation when considering how 
disk detectability evolves. 

Other loss processes are readily incorporated into this scheme. In Sj4] we also discussed 
loss rates that can be defined by a power law in particle size, including dynamical depletion, 
a more realistic prescription for radiation loss rates that are not simply inversely proportional 
to paticle size, and stellar wind drag. Steep power law loss rates lead to size distributions 
and evolutions that are directly analogous to P-R drag (but with modified exponents), but 
loss rates that are less biased to small particles are harder to treat. However, the application 
to dynamical depletion agrees with the results in the literature. The way in which P-R drag 
becomes independent of size for the smallest particles causes a ripple and flattening of this 
distribution at such sizes in the P-R drag regime. The application to stellar wind drag results 
in distributions and evolutions that are the same as those for P-R drag, but if this effect 
becomes important at all, it does so at higher disc masses meaning that it has the potential 
to be relevant for disks that are detectable with current technology. 
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There are now several surveys for debris disks that show how disk brightness drops with 
age (see review in Wyatt 2008), and further surveys are currently underway (e.g., Matthews 
et al. 2010). Interpretation of these surveys requires efficient schemes so that populations of 
large numbers of debris disks can be modelled (e.g., Wyatt et al. 2007b; Lohne et al. 2008). 
As such surveys push to lower disk masses, it is also important to be able to incorporate 
consideration of loss processes. The results presented here suggest that when disks drop to 
a level at which P-R drag becomes important, which is close to the limit of detectability 
with Spitzer (see Fig. 4 of Wyatt et al. 2007a), their brightness begins to drop rapidly, as 
fast as °c r~ 2 - 8 . Such a fast fall off is also expected for disks affected by stellar wind drag, 
which may occur for disks well above the current detectability threshold. In short, this paper 
provides an efficient model with which to interpret observations of the debris disks of nearby 
stars, and their evolution, including in regimes where loss processes are important for small 
(or large) particles. 
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